home *** CD-ROM | disk | FTP | other *** search
/ EnigmA Amiga Run 1996 February / EnigmA AMIGA RUN 04 (1996)(G.R. Edizioni)(IT)[!][issue 1996-02][Skylink CD III].iso / earcd / midi / gfft.lha / gfft-2.03 / source / gfft-2.03-source.lha / okspctrm.c < prev    next >
C/C++ Source or Header  |  1996-01-02  |  14KB  |  507 lines
  1. /***************************************************************************
  2.  *          Copyright (C) 1994  Charles P. Peterson                  *
  3.  *         4007 Enchanted Sun, San Antonio, Texas 78244-1254             *
  4.  *              Email: Charles_P_Peterson@fcircus.sat.tx.us                *
  5.  *                                                                         *
  6.  *          This is free software with NO WARRANTY.                  *
  7.  *          See gfft.c, or run program itself, for details.              *
  8.  *              Support is available for a fee.                      *
  9.  ***************************************************************************
  10.  *
  11.  * Program:     gfft--General FFT analysis
  12.  * File:        okspctrm.c
  13.  * Purpose:     OK! Do an fft spectrum analysis
  14.  * Author:      Charles Peterson (CPP)
  15.  * History:     24-July-1993 CPP; Created.
  16.  *               6-Feb-95 CPP (1.31); Progress Requester
  17.  *
  18.  * Comment:
  19.  *
  20.  *    As used here, a 'spectrum analysis' is one which starts with real
  21.  *    samples, and ends up with real magnitudes of some sort (possibly with
  22.  *    some sort of 'windowing' applied to the input which may be analyzed
  23.  *    in 'segments' which might be padded and/or overlapped, and possibly
  24.  *    with normalization applied to the output).  The output magnitudes may
  25.  *    either be a 'power' spectrum or an 'amplitude' or 'voltage' spectrum
  26.  *    (the latter being the positive square root of the former.)  *
  27.  *
  28.  *    (Ordinary fft, in contrast, may begin with complex input values and
  29.  *    usually returns complex output values.  Usually the input to an
  30.  *    ordinary fft is not segmented, and doing so may present additional
  31.  *    problems.  Most frequently, ordinary fft is actually part of some
  32.  *    larger problem, such as convolution or spectrum analysis.)
  33.  *
  34.  *    This function has gotten a little out of hand, primarily because of
  35.  *    all the possible combinations of options that are applied here.
  36.  *    Possibly 'Interleave' might be left out next time unless someone
  37.  *    notes a use for it.  (It didn't turn out to be as useful as I
  38.  *    expected.  When curves made with different amounts of interleave are
  39.  *    spliced together, they don't meet, even in 'PSDensity' mode.)  Likewise
  40.  *    for 'Pad,' an option which I've never seen do anything but make
  41.  *    matters much worse.  But somehow, some people seem to think that it's
  42.  *    useful, or at least harmless.  It isn't either, in my experience.
  43.  *    Leaving out those two options would reduce this size of this function
  44.  *    substantially.  
  45.  *
  46.  *    But, don't expect any such streamlining to make a signficant (or
  47.  *    even measureable) difference in performance.  My lprof tests indicate
  48.  *    that this entire beast takes less than 1% of the total time under
  49.  *    typical conditions.  Considering that there is real work going on
  50.  *    here--even if you don't use the more esoteric options--that is quite
  51.  *    reasonable.  (As expected, cfft takes around 70% of the time under
  52.  *    typical conditions, which isn't unreasonable either.)
  53.  *
  54.  *    Anyway, if you want to see my best programming and comment writing,
  55.  *    go to cfft.c, the function which actually does the fft, though it
  56.  *    has gotten harrier since I added an accelerated mode.
  57.  *
  58.  */
  59. #define INITIAL_READ_PROGRESS 5
  60. #define ENDING_WRITE_PROGRESS 45
  61.  
  62. #define USE_MEMMOVE 
  63. /* Also uses MEMCPY, where applicable */
  64. /* Measured with lprof, only about 0.3% improvement overall. */
  65. /* But, even with 1 bin, didn't make matters worse either. */
  66.  
  67. #include <time.h>  /* time and difftime */
  68. #include <string.h> /* memmove and memcpy */
  69.  
  70. #include "gfft.h"
  71. #include "settings.h"
  72. #include "wbench.h"    /* progress requester stuff */
  73.  
  74. double LoopTimeElapsed = 0.0;
  75.  
  76. double Percent_Per_Loop = 1.0;
  77. double Initial_Progress = (double) INITIAL_READ_PROGRESS;
  78. double Ending_Progress = (double) ENDING_WRITE_PROGRESS; /* Must not be 0 */
  79.  
  80. static BIN_TYPE *bins = NULL;
  81. static ULONG number_bins = 0;
  82. static int segment_count = 0;
  83.  
  84. void do_re_output (void)
  85. {
  86.     if (bins && number_bins && segment_count)
  87.     {
  88.     ok_write (bins, (long) number_bins, (long) segment_count);
  89.     }
  90.     else
  91.     {
  92.     error_message (CANT_RE_OUTPUT);
  93.     RAISE_ERROR (NOTHING_SPECIAL);    /* longjmp outa here */
  94.     }
  95. }
  96.  
  97.  
  98. ULONG ok_spectrum (BOOLEAN do_it_for_real)
  99. {
  100.     static float *indata = NULL;
  101.     static float *overlapdata = NULL;
  102.     static float *interleavedata = NULL;  /* e.g. Aseg Bseg Aseg Bseg ... */
  103.     ULONG actually_read = 0;
  104.     ULONG number_samples = 1;
  105.     ULONG number_interleave_samples = 1;
  106.     ULONG half_samples = 0;
  107.     int more_data = TRUE;
  108.     int error_in_loop = FALSE;
  109.     int overlap = Overlap;     /* May be changed if overlap impossible */
  110.     int last_partial_segment_used = FALSE;
  111.     ULONG loop_count = 0;
  112.     ULONG i;
  113.     ULONG j;
  114.     ULONG ai;
  115.     ULONG aj;
  116.     time_t time_beginning;
  117.     time_t time_ending;
  118.     ULONG total_actually_read;
  119.  
  120.     number_bins = 0;
  121.  
  122.     Sample_Sum_Of_Squares = 0.0;
  123.     Sample_Frame_Count = 0;
  124.  
  125.     total_actually_read = ok_read (NULL, NULL);
  126.  
  127.     if (NumberBins == MAX_BINS)
  128.     {
  129.     if (!total_actually_read)
  130.     {
  131.         error_message (NO_DATA_PRESENT);
  132.         RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
  133.     }
  134.     else
  135.     {
  136.         ULONG next_power;
  137.         ULONG total_in_leaf = total_actually_read / Interleave;
  138.         if (!total_in_leaf)
  139.         {
  140.         error_message (INSUFFICIENT_DATA);
  141.         RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
  142.         }
  143.         next_power = get_pos_power_2 (total_in_leaf);
  144.         if (Pad || next_power == total_in_leaf)
  145.         {
  146.         number_bins = next_power / 2; /* pad if required*/
  147.         }
  148.         else
  149.         {
  150.         number_bins = next_power / 4; /* truncate if required */
  151.         }
  152.         number_samples = number_bins * 2 * Interleave;
  153.         bins_d_message (total_actually_read, number_bins);
  154.     }
  155.     }
  156.     else
  157.     {
  158.     number_bins = (ULONG) NumberBins;
  159.     if (number_bins > 0)
  160.     {
  161.         number_samples = number_bins * 2 * Interleave;
  162.     }
  163.     else
  164.     {
  165.         number_samples = 1;
  166.     }
  167.     }
  168. /*
  169.  * Now we have number_samples and number_bins
  170.  */
  171.     if (!do_it_for_real) return number_bins;
  172.  
  173.     segment_count = 0;
  174.  
  175.     time_beginning = time (NULL);
  176.  
  177.     half_samples = number_samples / 2;
  178.     if (number_samples < 2)
  179.     {
  180. /*
  181.  * Can't overlap segments of length 1
  182.  */
  183.     overlap = FALSE;
  184.     }
  185.     if (bins)
  186.     {
  187.     gfree (bins);
  188.     }
  189.     if (indata)
  190.     {
  191.     gfree (indata);
  192.     }
  193.     if (overlapdata)
  194.     {
  195.     gfree (overlapdata);
  196.     overlapdata = NULL;
  197.     }
  198.     if (interleavedata)
  199.     {
  200.     gfree (interleavedata);
  201.     interleavedata = NULL;
  202.     }
  203.       
  204. /*
  205.  * Note that there is one extra bin for 0 Hz; bins must be pre-zeroed
  206.  */
  207.     bins = gcalloc ((number_bins + 1), (long) sizeof(BIN_TYPE), 
  208.             NOTHING_SPECIAL);
  209.     indata = gmalloc (number_samples * sizeof(float), 
  210.               NOTHING_SPECIAL);
  211.     if (overlap)
  212.     {
  213.     overlapdata = gmalloc (number_samples * sizeof(float), 
  214.                NOTHING_SPECIAL);
  215.     }
  216.     if (Interleave > 1)
  217.     {
  218.     number_interleave_samples = number_samples / Interleave;
  219.     if (!number_interleave_samples)
  220.     {
  221.         error_message (INSUFFICIENT_DATA);
  222.         RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
  223.     }
  224.     interleavedata = gmalloc (number_interleave_samples * sizeof(float),
  225.                  NOTHING_SPECIAL);
  226.     }
  227.  
  228. /*
  229.  * Now, figure the number of FFT Inner Loops (FILs) required
  230.  *   This is not going to deal with all possible complexities...
  231.  *   since there is no great harm if we are wrong--it's only for the
  232.  *   progress indicator.
  233.  *     In particular, it's not going to deal with Interleave and Pad,
  234.  *     which are fairly useless options not directly supported by the
  235.  *     GUI interface anyway.
  236.  */
  237.     {
  238.     extern long Actual_Time_Segments;
  239.     extern long Inner_Loop_Count;
  240.     long fils_per_segment;
  241.     long est_fils;
  242.     int segments_per_time_segment;
  243.     double computational_percentage;
  244.     int time_segments;
  245.  
  246.     fils_per_segment = fft_inner_loops (number_bins);
  247.  
  248.     if (!overlap)
  249.     {
  250.         segments_per_time_segment = total_actually_read /
  251.                                      number_samples;
  252.     }
  253.     else if (total_actually_read % number_samples == 0)
  254.     {
  255.         segments_per_time_segment = 2 * (total_actually_read / 
  256.                      number_samples) - 1;
  257.     }
  258.     else
  259.     {
  260.         segments_per_time_segment = 2 * (total_actually_read /
  261.                      number_samples);
  262.     }
  263.  
  264.     if (Time3D)
  265.     {
  266.         time_segments = Actual_Time_Segments;
  267.         Ending_Progress = 1;
  268.         Initial_Progress = 1;
  269.     }
  270.     else
  271.     {
  272.         Inner_Loop_Count = 0;
  273.         time_segments = 1;
  274.         if (segments_per_time_segment * time_segments > 2)
  275.         {
  276.         Ending_Progress = ENDING_WRITE_PROGRESS * 2 /
  277.           (segments_per_time_segment * time_segments);
  278.         Initial_Progress = INITIAL_READ_PROGRESS * 2 /
  279.           (segments_per_time_segment * time_segments);
  280.         if (Ending_Progress < 1) Ending_Progress = 1;
  281.         }
  282.         else
  283.         {
  284.         Ending_Progress = ENDING_WRITE_PROGRESS;
  285.         Initial_Progress = INITIAL_READ_PROGRESS;
  286.         }
  287.     }
  288.     
  289.     est_fils = fils_per_segment * segments_per_time_segment *
  290.                   time_segments;
  291.  
  292.     computational_percentage = ((100 - Initial_Progress) -
  293.                                 Ending_Progress);
  294.  
  295.     if (est_fils)
  296.     {
  297.         Percent_Per_Loop = computational_percentage / est_fils;
  298.     }
  299.     else
  300.     {
  301.         Percent_Per_Loop = computational_percentage;
  302.     }
  303.     if (number_bins > 1)
  304.     {
  305.         progress_requester__update (1);  /* Warm and fuzzy */
  306.     }
  307.     else
  308.     {
  309.         progress_requester__update (-1);  /* Can't really tell */
  310.     }
  311.  
  312.     }
  313. /*
  314.  * Now that these preliminaries are out of the way, we can
  315.  * begin looping over input data
  316.  */
  317.     time_beginning = time (NULL);
  318.  
  319.     while (more_data)
  320.     {
  321.     actually_read = ok_read (indata, number_samples);
  322.  
  323.       /* ignored if already > this value */
  324.     progress_requester__update ((int) Initial_Progress);
  325.  
  326.     if (!actually_read)
  327.     {
  328.         more_data = FALSE;
  329.         if (!loop_count)
  330.         {
  331.         error_message (NO_DATA_PRESENT);
  332.         error_in_loop = TRUE;
  333.         }
  334.         break;
  335.     }
  336.     if (actually_read < number_samples)  /* Deal with partial segment */
  337.     {
  338.         more_data = FALSE;
  339.         if (overlap && loop_count)
  340.         {
  341.         for (i = 0, j = actually_read; j < number_samples; 
  342.              i++, j++)
  343.         {
  344.             overlapdata[i] = overlapdata[j];
  345.         }
  346.         for (j = 0; j < actually_read; i++, j++)
  347.         {
  348.             overlapdata[i] = indata[j];
  349.         }
  350.         if (Interleave > 1)    /* Overlap and Interleave */
  351.         {
  352.             for (ai = 0; ai < Interleave; ai++)
  353.             {
  354.             float *a_data = interleavedata;
  355.             for (aj = ai; aj < number_samples; aj += Interleave)
  356.             {
  357.                 *a_data++ = overlapdata[aj];
  358.             }
  359.             ok_window (interleavedata, number_interleave_samples);
  360.             ok_rfft (interleavedata, number_interleave_samples);
  361.             ok_sigma (interleavedata, bins, number_bins);
  362.             segment_count++;
  363.             }
  364.         }
  365.         else  /* Overlap only...no Interleave */
  366.         {
  367.             ok_window (overlapdata, number_samples);
  368.             ok_rfft (overlapdata, number_samples);
  369.             ok_sigma (overlapdata, bins, number_bins);
  370.             segment_count++;
  371.         }
  372.         overlap = FALSE;  /* overlap buffer now used up */
  373.         last_partial_segment_used = TRUE;
  374.         }
  375.         if (Pad)
  376.         {
  377.         error_message (PADDING_TAILEND);
  378.         for (i = actually_read; i < number_samples; i++)
  379.         {
  380.             indata[i] = 0.0;
  381.         }
  382.         }
  383.         else
  384.         {
  385.         if (!segment_count)
  386.         {
  387.             error_message (INSUFFICIENT_DATA);
  388.             error_in_loop = TRUE;
  389.         }
  390. /*        else if (!last_partial_segment_used)
  391.         {
  392.             error_message (IGNORING_TAILEND);
  393.         } */
  394.         break;
  395.         }
  396.     }
  397. /*
  398.  * Here, we've got a full segment to work with (padded or not)
  399.  */
  400.     if (overlap)
  401.     {
  402.         if (loop_count)  /* Can't process overlap first time around */
  403.         {
  404. #ifdef USE_MEMMOVE
  405.         memmove (overlapdata, &overlapdata[half_samples], 
  406.              half_samples * sizeof (float));
  407.         memcpy (&overlapdata[half_samples], indata,
  408.             half_samples * sizeof (float));
  409. #else
  410.         for (i = 0, j = half_samples; i < half_samples; i++, j++)
  411.         {
  412.             overlapdata[i] = overlapdata[j];
  413.         }
  414.         for (j = 0; j < half_samples; i++, j++)
  415.         {
  416.             overlapdata[i] = indata[j];
  417.         }
  418. #endif
  419.         if (Interleave > 1)   /* Overlap and Interleave */
  420.         {
  421.             for (ai = 0; ai < Interleave; ai++)
  422.             {
  423.             float *a_data = interleavedata;
  424.             for (aj = ai; aj < number_samples; aj += Interleave)
  425.             {
  426.                 *a_data++ = overlapdata[aj];
  427.             }
  428.             ok_window (interleavedata, number_interleave_samples);
  429.             ok_rfft (interleavedata, number_interleave_samples);
  430.             ok_sigma (interleavedata, bins, number_bins);
  431.             segment_count++;
  432.             }
  433.         }
  434.         else  /* Overlap only...no Interleave */
  435.         {
  436.  
  437.             ok_window (overlapdata, number_samples);
  438.             ok_rfft (overlapdata, number_samples);
  439.             ok_sigma (overlapdata, bins, number_bins);
  440.             segment_count++;
  441.         }
  442.         }
  443. #ifdef USE_MEMMOVE
  444.         memcpy (overlapdata, indata, number_samples * sizeof (float));
  445. #else
  446.         for (i = 0; i < number_samples; i++)
  447.         {
  448.         overlapdata[i] = indata[i];
  449.         }
  450. #endif
  451.     }
  452.     if (Interleave > 1)  /* Regular (unoverlapped) alternation */
  453.     {
  454.         for (ai = 0; ai < Interleave; ai++)
  455.         {
  456.         float *a_data = interleavedata;
  457.         for (aj = ai; aj < number_samples; aj += Interleave)
  458.         {
  459.             *a_data++ = indata[aj];
  460.         }
  461.         ok_window (interleavedata, number_interleave_samples);
  462.         ok_rfft (interleavedata, number_interleave_samples);
  463.         ok_sigma (interleavedata, bins, number_bins);
  464.         segment_count++;
  465.         }
  466.     }
  467.     else
  468.     {
  469.         ok_window (indata, number_samples);
  470.         ok_rfft (indata, number_samples);
  471.         ok_sigma (indata, bins, number_bins);
  472.         segment_count++;
  473.     }
  474.     loop_count++;
  475.     }
  476.     time_ending = time (NULL);
  477.     LoopTimeElapsed = difftime (time_ending, time_beginning);
  478.  
  479.     if (!error_in_loop)
  480.     {
  481.     loop_time_message (LoopTimeElapsed);
  482.     ok_write (bins, number_bins, segment_count);
  483.     }
  484.     gfree (indata);
  485.     indata = NULL;
  486. /*    gfree (bins); save for ReOut */
  487. /*    bins = NULL; */
  488.     if (overlapdata)
  489.     {
  490.     gfree (overlapdata);
  491.     overlapdata = NULL;
  492.     }
  493.     if (interleavedata)
  494.     {
  495.     gfree (interleavedata);
  496.     interleavedata = NULL;
  497.     }
  498.     if (error_in_loop)
  499.     {
  500.     RAISE_ERROR (NOTHING_SPECIAL);  /* longjmp outa here! */
  501.     }
  502.  
  503.     return number_bins;
  504. }
  505.  
  506.     
  507.